Comparative simulation study of colloidal gels and glasses 
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Using computer simulations, we identify the mechanisms causing aggregation and structural 
arrest of colloidal suspensions interacting with a short-ranged attraction at moderate and high 
densities. Two different non-ergodicity transitions are observed. As the density is increased, a glass 
transition takes place, driven by excluded volume effects. In contrast, at moderate densities, gelation 
is approached as the strength of the attraction increases. At high density and interaction strength, 
both transitions merge, and a logarithmic decay in the correlation function is observed. All of these 
features are correctly predicted by mode coupling theory. 



Colloidal dispersions aggregate into various non- 
equilibrium structures depending on density, interaction 
strength and range. The accompanying rheology and 
structure are among the key properties desired for their 
technological applications [|]. Moreover, thanks to the 
possibility to tailor effective interactions by e.g. addition 
of salt and polymer, colloids allow us to study the funda- 
mental mechanisms of kinetic arrest. Whereas colloidal 
hard spheres have become a model system for the study 
of structural arrest at a glass transition Q , colloidal gela- 
tion has only recently been associated with glassy behav- 
ior 1^-^. Colloidal gelation is ubiquitous in suspensions 
driven by attractions of quite short range and moderate 
to high strength ||^. At low packing fractions, it entails 
the formation of heterogeneous and often self-similar net- 
works; there, an interplay of phase separation kinetics 
and percolation often are considered responsible for its 
existence At higher densities, the gelation boundary 
extends into the homogeneous fluid region §,0, where 
it also lies well separated from estimates of percolation 
S0~P ■ Crossing into the gelled state anywhere along the 
transition line results in qualitatively the same phenom- 
ena, like flow properties that indicate solidification |0, |j , 
and non-ergodic dynamics according to light scattering 

We present simulations designed to identify the mech- 
anism of colloidal gelation driven by attractions of only 
moderate strength. Because of the distance of the gel 
boundary from other boundaries (percolation and phase 
separation) at higher densities, we concentrate on these, 
where we sweep out the region between gel and glass 
transition lines. We show that both non-equilibrium 
transitions are caused by a slowing down of local re- 
arrangements, as predicted well by mode coupling the- 
ory (MCT) ||,|l^,|l3|. We contrast the glass transition, 
caused by caging of particles owing to steric hindrance, 
with attraction-driven gelation caused by bonding be- 
tween particles. We verify that the simultaneous pres- 
ence of two non-ergodic states results in anomalous non- 
exponential (logarithmic) time dependences, as recently 



conjectured to explain observations in micellar systems 
1^ or microgel suspensions (ij] . 

The simulated system comprises 1000 soft-core {V{r — 
\ri — Yj\) oc {uij/r)^^, ttij = ai + aj) particles of mean ra- 
dius a with polydispersity in size (flat distribution with 
10% width) to prevent crystallisation. Densities are re- 



ported as packing fractions 



-na . A short range 



attraction, mimicking the polymer induced depletion at- 
traction in experimental systems j^,^,|l^, is given by 
an Asakura-Oosawa (AO) form generalized to polydis- 
perse systems [^,|l^. The range of the attraction, 2^, is 
set to 0.2a, and its strength is proportional to the poly- 
mer concentration (pp. To help avoid liquid-gas separa- 
tion, a weak long range barrier is added to the potential. 
The barrier extends from ai2 + 2^ to 4a, and is described 
by a fourth order polynomial matched to give a continu- 
ous force. Its maximal height is IkgT, which equals the 



depth at contact of the AO potential at 



0.0625. In 



all states studied, the barrier is much smaller than the 
attraction, and in the purely repulsive case {4>p — 0) it 
is omitted. We will measure lengths in units of a and 
time in units of \J jZv^ , where the thermal velocity, 
V, is set to 2/\/2,. Equations of motion were integrated 
using the velocity- Verlet algorithm, with a time step of 
0.0025. Colloidal dynamics (neglecting hydrodynamic in- 
teractions) were mimicked by running the simulations in 
the canonical (constant NTV) ensemble, where the ther- 
mostat plays the role of the surrounding liquid. Every 
N time steps, the velocity of the particles was rescaled 
to assure constant v. No effect of TV on the results was 
observed for well equilibrated samples. 

The central quantity of our study will be the self 
part of the intermediate scattering function, $^(i) = 
(expiq-(rj(t) — rj(0))), for wave- vector q, where (. . .) de- 
notes an average over particles and time origin, '^^{t) al- 
lows us to probe and identify the nature of the dominant 
dynamical mechanism because of (i) its g-dependence 
and (m) the detailed predictions that are available from 
MCT. Indeed, if a structural arrest at a non-ergodicity 
transition is approached, <i>^(i) reveals a two-step pro- 
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cess, where the decay from the plateau is given by the 
von Schweidler power-law series p7|] : 



Here /| is the non-ergodicity parameter, h^p and ftg^^ 
are amplitudes, and 6 is known as the von Schweidler 
exponent. On the one hand, the observation of this (uni- 
versal) von Schweidler behavior - and tests of further 
relations, as done below - establishes that a feed-back 
mechanism in the structural relaxation causes arrest. On 
the other hand, the (non-universal) wave-vector depen- 
dence of the amplitudes, like /|, allows us to identify the 
specific kinetic process which freezes out. As the transi- 
tion is approached, the characteristic time r diverges as 
T <x.\4> — (t>'^\^^ where 7 is determined by 6, see e.g. [nTl. 
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FIG. 1. Relaxation time r as a function of 



0.45 
for three 



colloid volume fractions: ipc 
4>c = 0.55 (O)- Inset: r vs. 



0.40(0), <j,c = 0.50(D), 
for soft spheres {4>p = 0: 



V{r) 



Dotted lines are fittings with predetermined 7. 



Figure 1 presents evidence for both the repulsion and 
attraction driven glass transitions, as identified by a di- 
verging T jl^ . Upon increasing the packing fraction ip^ 
(inset in figure 1), the system approaches a glass tran- 
sition caused by steric hindrance, which we have stud- 
ied including only the r^^^ repulsion (0p = and no 
barrier) jl^. The transition correlates well with obser- 
vations at the colloidal glass transition j2| and previous 
simulations of e. g. a glassy Lennard- Jones mixture [ pO[ . 
We have analysed it using the concepts of idealized MCT, 
but will present only a few results here for comparison 
with gelation. Gelation itself is induced by strengthening 
the attraction at intermediate packing fractions (figure 
1). There, far from equilibrium or percolation transi- 
tions (as we tested by monitoring the static structure 
factor), arrest again is of kinetic origin, and occurs at 
lower attraction strengths the higher ^c- 
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FIG. 2. (a): Correlation functions at (pc = 0.4, 
(j)p — 0.425 and von Schweidler fits. From top to botton, 
q = 3.9, 6.9, 9.9, 15, 20, 25, 30, 35, 40, 50. (b): /| for glass 
(o) and gel (•) (from (a)) transitions, with the Gaussian ap- 
proximation for both of them (dotted lines) and the MCT 
result for hard spheres (dashed line) | 



Inset: /ig^' using r 



from figure 1, and MCT result for hard spheres. 

To shed light on this transition, the correlation func- 
tions at different wave-vectors were studied. The slowest 
state at cjic — 0.40 is presented in Fig |[ The upper panel 
shows the self intermediate scattering functions for dif- 
ferent wave-vectors, and the fits using (|l|) up to second 
order. A common exponent b was taken in the fitting, 
yielding a value of 6 = 0.38, appreciably lower than the 
hard spheres value b — 0.53 (which we found for our soft 
sphere glass at (j)p — 0). As predicted by MCT, we can 
calculate from b the divergence of the relaxation times in 
Fig. |l|. The resulting value 7 — 3.03 fits the data, while 
7 = 2.63 at the soft sphere glass transition The gel 
transition is estimated to occur at 4>p — 0.431. Because 
we find the universal properties predicted by MCT, we 
conclude that at (j>c — 0.40 the gel transition is a regular 
non-ergodicity transition in the structural dynamics. 

Their very different g-width (Fig. ^) for the non- 
ergodicity parameters and amplitudes brings out a major 
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difference in the two underlying mechanisms. Whereas 
repulsions localize the particle within a cage, which it 
can explore up to mean squared displacements rf of the 
order of rf = 0.13 (from our simulations, not shown; 
rf = 0.134 from MCT attractions bind the particle 

to its neighbors and thus localize it much more tightly. 
At (/)c — 0.40, we find rf = 0.018 by simulations, which 
is of the order of a low-density estimate [Q for our in- 
teraction range, 2^ ~ 0.2. The corresponding high am- 
plitudes /| of density fluctuations are consistent with 
light scattering observations at fixed q |l^, ^ and with 
MCT calculations 0,|l|,|l^. These fluctuations extend 
to large q and relax only when the particles break free 
from their bonds. The comparison with the Gaussian 
approximation, /^^ = exp {— g^r;^/6} evidences stronger 
non-Gaussian effects at gelation than at the glass tran- 
sition. We stress the cooperativity of the structural re- 
laxation at both transitions. Holding all particles fixed, 
except for one, leads to mean squared displacement for 
the tracer (as it explores the frozen enviroment) much 
smaller than in the free system (before the start of the 
structure relaxation of the free system, the ratio is ~ 6 
for both cases) . The cage or network of bonds around an 
arrested particle thus necessarily fluctuates with it. 




FIG. 3. Correlation functions vs. rescaled time t/r for 
(pc = 0.4 and (pp = 0.375, 0.39, 0.40, 0.415, 0.425 (from right 
to left) and von Schweidler fit. Inset: Same plot for soft 
spheres and (j>c = 0.53,0.54,0.55,0.56,0.57, 0.58,0.585,0.59, 
and MCT master curve Q. 

To test further the nature of the gel transition, the 
scaling of the final (or a-) decay was studied. In Fig. ^ 
we present the rescaled ($^(i/T = 1) = 0.25) correlation 
functions at q — 9.9 for different attraction strengths, 
close to the gel transition. In the inset to this figure, a 
similar plot deals with the glass transition {q — 3.9 in this 
case). In both cases, the curves clearly collapse during 
the a-decay indicating an unique mechanism which dom- 
inates the slowing down at the transitions. For the purely 
repulsive case, the MCT master curve for the rescaled de- 



cay of hard spheres at a slightly larger wave-vector 
(g = 4.3) is also presented, confirming the quantitative 
agreement between MCT and our results. In the gel case, 
no master function is available, but the fit to (|^) is pre- 
sented. The different stretching in the two cases is clear. 
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FIG. 4. Same as figure 2(a) for (j>c ~ 0.55 and (jjp = 0.375. 
Dotted lines: logarithmic fittings to the correlators around 
fq. The dashed line represents ti. 

We study now the gel transition at a higher density 
(j)c = 0.55, where it lies closer to the glass transition. 
Within MCT, the simultaneous existence of two differ- 
ent non-ergodicity transition branches opens the possibil- 
ity for end-point singularities, where the branches merge 
|2l| . In systems with short range attractions, such a sin- 
gularity has been predicted close to the crossing of the 
two transition lines (the actual distance depending on 
the details of the potential) |lj, ^ . Close to the sin- 
gularity, a logarithmic decay around the plateau in the 
correlation function, ^^{t) = f^"^ — Cq\n[t/ti), is a pro- 
posed signature ||l^, Q , and intriguingly is observed ex- 
perimentally in more complicated systems Q . Having 
identified two different non-ergodicity transitions in our 
system, we now test this prediction. In Fig. ^ the corre- 
lation functions at the same wave- vectors as in Fig. |^ are 
presented for the state 0p = 0.375. Logarithmic decays 
are observed in all of the correlators (linear traces in the 
plot) for up to three decades in time, signalling a higher 
order singularity nearby. It is interesting to note that 
the logarithmic trace of the correlator has different ex- 
tents, depending on the wave-vector. To make it clearer, 
{'^q ~ fq^)/C^q ^s. timc has been plotted as an inset to 
Fig. H where fq^ is determined at ii = 5 (vertical line 
in Fig. ^). Deviations from the logarithmic decay are 
stronger, the higher /|. This is in complete agreement 
with the theoretical expectations in |l3| . Since both the 
von Schweidler decay (associated with the gel transition) 
and the logarithmic trace now take place in the same 
window, important corrections to the a-scaling of the 
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curves are expected. This is seen in the main graph of 
Fig. ||: the long time dynamics at different states can- 
not be collapsed onto a master curve by time rescaling. 
This shows that two mechanisms are responsible for the 
structural slowing down and that changes in the control 
parameters change the relative distance to gelation but 
also to the higher order singularity. The intermediate 
isochore, (j)c = 0.50, shows a mixed behavior: a logarith- 
mic decay in a smaller window than for (pc — 0.55, and 
followed by an apparent power law decay. 
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FIG. 5. Correlation function vs. 
q = 20. (t>c = 0.55, and 4)p = 0.325, 



rescaled time, t/r, at 
0.34, 0.35, 0.365, 0.375 



from right to left at t < r. Inset: ($g — fq)/Cq as a func- 
tion of time for the same state and wave- vectors as figure |^ 
(wave- vectors increasing from top to botton). 

In summary, using MD simulations, we have deduced 
from the wave-vector dependence of the dynamical den- 
sity fluctuations, that repulsion and short-ranged at- 
traction lead to two different structural arrests at high 
enough density or attraction strength, respectively. At 
the merging of both glassy states, subtle logarithmic time 
variations appear. Comparing with the recent MCT pre- 
dictions of these phenomena we find perfect agreement. 
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